Introdução a Stan

(e talvez RSample)

Carolina Musso

Estatística Computacional

Roteiro da apresentação

  • Inferência bayesiana numa casca de noz
    • Conceitos centrais, vantagens e problemas
  • Introdução STAN
    • Programação probabística
  • Como usar o STAN no R
    • Rstan, rstanarm, tidymodels
    • explorar os resultados com ShinySTAN
  • RSample
    • num sei se vai chegar até aqui…

Bayesiana x Frequentistas

Uma outra forma de enxergar a probabilidade/incerteza

Frequentistas

  • Parâmetro fixo.
  • Amostragens repetidas teoricamente
  • Intervalo de confiança: 95% dos intervalos irão conter o parâmetro.

Bayesianos

  • Incerteza sobre o parâmetro
  • v.a. da perspectiva do observador
  • Intervalo de credibilidade: Há 95% de chance do parâmetro estar nesse intervalo.
  • ZedStatistics - Canal bem nteressante com interpretações visuais e intuitivas de estatística.

Teorema de Bayes

Bayes/Laplace: séc. 18 e 19

  • \[ P (A|B) = \frac{P(B|A)*P(A)}{P(B)} \]

  • P(B|A): verossimilhança

  • P(A): incerteza sobre o parâmetro (priori)

  • P(B) é a constatante normalizadora

Resumindo o raciocínio

\[ \pi(\theta|y) \propto \pi(y|\theta) * \pi(\theta) \]

\[ posteriori = verossimilhança * priori \]

  • Consigo atualizar minhas estimativas e interpretar os resultados de forma mais direta

    • Não tem almoço gratis

      • Definir a priori é subjetivo.

      • Computar a posteriori as vezes é bem difícil

      • Maldição da dimensionalidade

Visualizando essa relação

Um exemplo inédito: lançamento de uma moeda

  • Acho balenceada, mas em 10 lançamentos teve 7 caras.

Para fazer modelagens bayesianas

  • Escrever um script que defina um modelo bayesiano

  • Conseguir amostrar e medir parâmetros da distribuição posteriori

    • STAN, 2012 (atualmente v.2.6)

      • Stanislaw Ulam, matetático polonês que desenvolveu o método de Monte Carlo nos anos 40.

      • “Sampling Through Adaptive Neighborhoods”

STAN

  • Software em C++ gratuito de código aberto

  • Capaz de

    • Inferência Bayesiana com amostrador MCMC (NUTS, HMC)

    • Inferência Bayesiana aproximada com “inferência variacional” (ADVI)

    • Estimativa de máxima verossimilhança penalizada com otimização (L-BFGS)

Programação probabilística

  • maioria das aplicações estatísticas.

  • Pode-se implementar automaticamente condicionamento.

STAN

NUTS

Lógica do STAN

\[ \pi(y,\theta) = \pi(y|\theta) * \pi(\theta) \]

  • Encontrar a distribuição conjunta

  • A linguagem de modelagem especifica os elementos em blocos:

    • dados

    • parâmetros

    • modelo

  • Palestra bem completa do Michael Betancour

Lógica do STAN

Primeiro bloco: espeço observado

data {                      // Data block
  int<lower = 1 >  N;       // Sample size
  int<lower = 1> K;         // Dimension of model matrix
  matrix[N, K] X;           // Model Matrix
  vector[N] y;              // Target variable
}

Segundo bloco: parâmetros

parameters {                // Parameters block}
  vector[K] beta;           // Coefficient vector
  real<lower = 0> sigma;    // Error scale, lower bound at 0
}

Lógica do STAN

Terceiro bloco: modelo

model {  
vector[N] mu;
  mu = X * beta;            // Creation of linear predictor
  
  // priors
  beta  ~ normal(0, 10);  
  sigma ~ cauchy(0, 5);
  
  // likelihood
  y ~ normal(mu, sigma);
  }

Tudo junto

data {                      
  int<lower = 1 >  N;       
  int<lower = 1> K;        
  matrix[N, K] X;           
  vector[N] y;              
}
parameters {                
  vector[K] beta;           
  real<lower = 0> sigma;    
}
model {                     
  vector[N] mu;
  mu = X * beta;            
  beta  ~ normal(0, 10);  
  sigma ~ cauchy(0, 5);
  y ~ normal(mu, sigma);
}
  • um arquivo .stan ou como um string

Código é compilado

  • Declarar as variáveis

  • Os valores específicos serão fornecidos pelos algortimos que avaliam o Stan.

  • O codigo é traduzido pra C++ e compilado.

    • maior eficiência, sem ter que programar em C++

  • Bibliotecas convertem em um programa executável capaz de avaliar a log-posteriori e a função gradiente

RSTan: instalar C++ toolchain

R versão >= 3.4.0 or later, de preferência >= 4.0.0

Recomenda-se RStudio version >= 1.4.

Necessário configurar o R para ser capaz de compilar código C++ code: Windows, Mac, Linux.

install.packages("rstan", repos = "https://cloud.r-project.org/", 
                 dependencies = TRUE)

library("rstan")
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

Exemplos

Mão na massa

Simular alguns dados …

set.seed(42)
N <- 1000; 
K = 1  # se quiser colocar mais covariaveis
coefs = c(5, .2); sigma = 2

covariates = replicate(
  K, 
  rnorm(n = N)) # se quiser colocar mais covariaveis

x <- rnorm(n = N); X = cbind(Intercept = 1,covariates)

mu = X %*% coefs
y <- rnorm(N, mu, sigma)

Usando lm()

modlm = lm(y ~ ., data = data.frame(X[, -1]))
summary(modlm)

Call:
lm(formula = y ~ ., data = data.frame(X[, -1]))

Residuals:
    Min      1Q  Median      3Q     Max 
-6.3624 -1.2669  0.0346  1.3258  6.9920 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.99245    0.06516  76.624   <2e-16 ***
X....1.      0.15604    0.06500   2.401   0.0166 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.06 on 998 degrees of freedom
Multiple R-squared:  0.005741,  Adjusted R-squared:  0.004745 
F-statistic: 5.763 on 1 and 998 DF,  p-value: 0.01655

Usando RStan

Dados devem ser uma lista!

X = cbind(Intercept = 1,x)
dados_stan = list(
  N = N,
  K = ncol(X),
  y = y,
  X = X
)

# outra opção
dados_stan <- read_rdump("dados.R")

Ajustando o modelo

fit_Rstan <- stan(
  model_code = stanmodelcode, # é um objeto no R criado com uma string 
  data   = dados_stan
)

# OU

fit_Rstan <- stan(
  file='fit_data.stan', # o programa stan está em um arquivo separado
  data   = dados_stan)

Vai demorar um pouquinho …


SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 7.7e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.77 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 0.089 seconds (Warm-up)
Chain 1:                0.086 seconds (Sampling)
Chain 1:                0.175 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 4.5e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.45 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 0.086 seconds (Warm-up)
Chain 2:                0.089 seconds (Sampling)
Chain 2:                0.175 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 4.6e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.46 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 0.087 seconds (Warm-up)
Chain 3:                0.092 seconds (Sampling)
Chain 3:                0.179 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 2.4e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.24 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 0.091 seconds (Warm-up)
Chain 4:                0.089 seconds (Sampling)
Chain 4:                0.18 seconds (Total)
Chain 4:

RStan

print(
  fit_Rstan,
  pars   = c('beta', 'sigma'),
  digits = 3,
  prob   = c(.025, .5, .975)
)
Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

         mean se_mean    sd  2.5%   50% 97.5% n_eff  Rhat
beta[1] 4.991   0.001 0.064 4.866 4.992 5.118  3651 0.999
beta[2] 0.218   0.001 0.063 0.091 0.218 0.336  3870 1.000
sigma   1.975   0.001 0.044 1.894 1.974 2.062  3880 1.001

Samples were drawn using NUTS(diag_e) at Sat Feb 11 22:44:02 2023.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

Alguns parametros interessantes

stan(..., 
  pars = NA, #se quiser salavr só alguns parâmetros
  chains = 4, # numero de cadeiras paralelsas
  iter = 2000, # numero de iterações
  warmup = floor(iter/2), # burnin
  thin = 1, # thining (período para se salvar amostras)
  init = "random", # valores iniciais
  seed = sample.int(.Machine$integer.max, 1), 
  algorithm = c("NUTS", "HMC", "Fixed_param"),
  ...)

Usando rstanarm

Uma interface da interface?

pacman::p_load(rstanarm)

beta_prior <- rstanarm::normal(0,10); sigma_prior <- rstanarm::cauchy(0,5)

fit_rstanarm<-rstanarm::stan_glm(y~., data=data.frame(X[, -1]), 
                         prior=beta_prior,
                         prior_aux = sigma_prior)

Coeficientes

fit_rstanarm$coefficients
(Intercept)     X....1. 
  4.9879549   0.2202713

Tidymodels

“The tidymodels framework is a collection of packages for modeling and machine learning using tidyverse principles”

Com o tidymodels

library(tidymodels)

bayes_mod <-   
  linear_reg() %>% 
  set_engine("stan", 
             prior_intercept = beta_prior, 
             prior = beta_prior,
             prior_aux = sigma_prior) 

bayes_fit <- 
  bayes_mod %>% 
  fit(y ~ x, data = data.frame(X[, -1]))

print(bayes_fit, digits = 5)
parsnip model object

stan_glm
 family:       gaussian [identity]
 formula:      y ~ x
 observations: 1000
 predictors:   2
------
            Median  MAD_SD 
(Intercept) 4.99110 0.06260
x           0.21982 0.05913

Auxiliary parameter(s):
      Median  MAD_SD 
sigma 1.97422 0.04356

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

Shinystan

pacman::p_load(shinystan)

avaliando_modelo <- launch_shinystan(fit_Rstan)

Outras referências

Um estudo de caso, Rbloggers, Outro, Stan for epdemiology

Become a Baysian with R and Stan, Bayesian data analysis - RStan demos

R/Stan examples, Stan function references

Prior distributions for RStan models

Materiais da StanCon, Lista de Estudos de Casos - oficial

Towards A Principled Bayesian Workflow

Bayesian Modeling Using Stan

Um dos pacotes do tidymodels é o RSample

O pacote rsample fornece funções para criar diferentes tipos de reamostragens com o intuito de:

  • reamostrar para estimar a distribuição de uma estatística.
  • estimar a perfomance de um modelo usando um conjunto de validação

Rsample

Os datasets criados são diretamente acessíveis no objeto mas sem consumir tanta memória.

library(rsample)
library(mlbench)

data(LetterRecognition)
lobstr::obj_size(LetterRecognition)
2.64 MB
set.seed(35222)
boots <- bootstraps(LetterRecognition, times = 50)
lobstr::obj_size(boots)
6.69 MB
# Object size per resample
lobstr::obj_size(boots)/nrow(boots)
133.74 kB
# Fold increase is <<< 50
as.numeric(lobstr::obj_size(boots)/lobstr::obj_size(LetterRecognition))
[1] 2.528489

Rsample

Objeto da classe rset`

library(rsample)
set.seed(8584)
bt_resamples <- bootstraps(mtcars, times = 3)
bt_resamples
# Bootstrap sampling 
# A tibble: 3 × 2
  splits          id        
  <list>          <chr>     
1 <split [32/14]> Bootstrap1
2 <split [32/12]> Bootstrap2
3 <split [32/14]> Bootstrap3

Cada amostra

first_resample <- bt_resamples$splits[[1]]
first_resample
<Analysis/Assess/Total>
<32/14/32>
#> <Analysis/Assess/Total>
#> <32/14/32>

as.data.frame(first_resample)
                     mpg cyl  disp  hp drat    wt  qsec vs am gear carb
Fiat 128            32.4   4  78.7  66 4.08 2.200 19.47  1  1    4    1
Toyota Corolla      33.9   4  71.1  65 4.22 1.835 19.90  1  1    4    1
Toyota Corolla.1    33.9   4  71.1  65 4.22 1.835 19.90  1  1    4    1
AMC Javelin         15.2   8 304.0 150 3.15 3.435 17.30  0  0    3    2
Valiant             18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
Merc 450SLC         15.2   8 275.8 180 3.07 3.780 18.00  0  0    3    3
Porsche 914-2       26.0   4 120.3  91 4.43 2.140 16.70  0  1    5    2
AMC Javelin.1       15.2   8 304.0 150 3.15 3.435 17.30  0  0    3    2
Dodge Challenger    15.5   8 318.0 150 2.76 3.520 16.87  0  0    3    2
Merc 450SLC.1       15.2   8 275.8 180 3.07 3.780 18.00  0  0    3    3
Porsche 914-2.1     26.0   4 120.3  91 4.43 2.140 16.70  0  1    5    2
Maserati Bora       15.0   8 301.0 335 3.54 3.570 14.60  0  1    5    8
Toyota Corona       21.5   4 120.1  97 3.70 2.465 20.01  1  0    3    1
Valiant.1           18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
Porsche 914-2.2     26.0   4 120.3  91 4.43 2.140 16.70  0  1    5    2
Lincoln Continental 10.4   8 460.0 215 3.00 5.424 17.82  0  0    3    4
Dodge Challenger.1  15.5   8 318.0 150 2.76 3.520 16.87  0  0    3    2
Mazda RX4           21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
Hornet Sportabout   18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2
Camaro Z28          13.3   8 350.0 245 3.73 3.840 15.41  0  0    3    4
Datsun 710          22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
Pontiac Firebird    19.2   8 400.0 175 3.08 3.845 17.05  0  0    3    2
Fiat 128.1          32.4   4  78.7  66 4.08 2.200 19.47  1  1    4    1
Ferrari Dino        19.7   6 145.0 175 3.62 2.770 15.50  0  1    5    6
Valiant.2           18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
Toyota Corolla.2    33.9   4  71.1  65 4.22 1.835 19.90  1  1    4    1
Pontiac Firebird.1  19.2   8 400.0 175 3.08 3.845 17.05  0  0    3    2
Duster 360          14.3   8 360.0 245 3.21 3.570 15.84  0  0    3    4
Mazda RX4.1         21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
Ford Pantera L      15.8   8 351.0 264 4.22 3.170 14.50  0  1    5    4
Duster 360.1        14.3   8 360.0 245 3.21 3.570 15.84  0  0    3    4
AMC Javelin.2       15.2   8 304.0 150 3.15 3.435 17.30  0  0    3    2

Analysis e Assessment

head(analysis(first_resample))
                  mpg cyl  disp  hp drat    wt  qsec vs am gear carb
Fiat 128         32.4   4  78.7  66 4.08 2.200 19.47  1  1    4    1
Toyota Corolla   33.9   4  71.1  65 4.22 1.835 19.90  1  1    4    1
Toyota Corolla.1 33.9   4  71.1  65 4.22 1.835 19.90  1  1    4    1
AMC Javelin      15.2   8 304.0 150 3.15 3.435 17.30  0  0    3    2
Valiant          18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
Merc 450SLC      15.2   8 275.8 180 3.07 3.780 18.00  0  0    3    3
head(assessment(first_resample))
                mpg cyl  disp  hp drat    wt  qsec vs am gear carb
Mazda RX4 Wag  21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
Hornet 4 Drive 21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
Merc 240D      24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2
Merc 230       22.8   4 140.8  95 3.92 3.150 22.90  1  0    4    2
Merc 280       19.2   6 167.6 123 3.92 3.440 18.30  1  0    4    4
Merc 280C      17.8   6 167.6 123 3.92 3.440 18.90  1  0    4    4

Intervalo de confiança com Bootstrap

pacman::p_load(nlstools,GGally)
data(O2K)

ggplot(O2K, aes(x = t, y = VO2)) + 
  geom_point()+
  theme_bw(base_size = 14)

O modelo

nonlin_form <-  
  as.formula(
    VO2 ~ (t <= 5.883) * VO2rest + 
      (t > 5.883) * 
      (VO2rest + (VO2peak - VO2rest) * (1 - exp(-(t - 5.883) / mu)))
    )

# valores iniciais por inspeção visual
start_vals <- list(VO2rest = 400, VO2peak = 1600, mu = 1)

res <- nls(nonlin_form, start = start_vals, data = O2K) 

tidy(res)
# A tibble: 3 × 5
  term    estimate std.error statistic  p.value
  <chr>      <dbl>     <dbl>     <dbl>    <dbl>
1 VO2rest   357.     11.4         31.3 4.27e-26
2 VO2peak  1631.     21.5         75.9 1.29e-38
3 mu          1.19    0.0766      15.5 1.08e-16

Usando cada amostra

# função para ajutar cada amostra ao modelo 
fit_fun <- function(split, ...) {
  nls(nonlin_form, data = analysis(split), ...) %>%
    tidy()
}

set.seed(462)
nlin_bt <-
  bootstraps(O2K, times = 2000, apparent = TRUE) %>%
  mutate(models = map(splits, ~ fit_fun(.x, start = start_vals)))
nlin_bt
# Bootstrap sampling with apparent sample 
# A tibble: 2,001 × 3
   splits          id            models          
   <list>          <chr>         <list>          
 1 <split [36/14]> Bootstrap0001 <tibble [3 × 5]>
 2 <split [36/13]> Bootstrap0002 <tibble [3 × 5]>
 3 <split [36/16]> Bootstrap0003 <tibble [3 × 5]>
 4 <split [36/12]> Bootstrap0004 <tibble [3 × 5]>
 5 <split [36/16]> Bootstrap0005 <tibble [3 × 5]>
 6 <split [36/13]> Bootstrap0006 <tibble [3 × 5]>
 7 <split [36/15]> Bootstrap0007 <tibble [3 × 5]>
 8 <split [36/16]> Bootstrap0008 <tibble [3 × 5]>
 9 <split [36/11]> Bootstrap0009 <tibble [3 × 5]>
10 <split [36/13]> Bootstrap0010 <tibble [3 × 5]>
# … with 1,991 more rows
nlin_bt$models[[1]]
# A tibble: 3 × 5
  term    estimate std.error statistic  p.value
  <chr>      <dbl>     <dbl>     <dbl>    <dbl>
1 VO2rest   359.      10.7        33.5 4.59e-27
2 VO2peak  1656.      31.1        53.3 1.39e-33
3 mu          1.23     0.113      10.9 2.01e-12

Obserservando os resultados

nls_coef <- 
  nlin_bt %>%
  dplyr::select(-splits) %>%
  # empilhar as colunas 
  unnest() %>%
  dplyr::select(id, term, estimate) 
head(nls_coef)
# A tibble: 6 × 3
  id            term    estimate
  <chr>         <chr>      <dbl>
1 Bootstrap0001 VO2rest   359.  
2 Bootstrap0001 VO2peak  1656.  
3 Bootstrap0001 mu          1.23
4 Bootstrap0002 VO2rest   358.  
5 Bootstrap0002 VO2peak  1662.  
6 Bootstrap0002 mu          1.26

Só os histogramas

p_ints <- int_pctl(nlin_bt, models)
nls_coef %>% 
  ggplot(aes(x = estimate)) + 
  geom_histogram(bins = 20, col = "white") + 
  facet_wrap(~ term, scales = "free_x") + 
  geom_vline(data = p_ints, aes(xintercept = .lower), col = "red") + 
  geom_vline(data = p_ints, aes(xintercept = .upper), col = "red")

Obrigada!

só de pensar que ainda tem uma lista pra entregar…